home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / frasrc19.zip / FRACTALB.C < prev    next >
C/C++ Source or Header  |  1995-03-04  |  27KB  |  991 lines
  1. /* -----------------------------------------------------------------
  2.  
  3. This file contains the "big number" high precision versions of the
  4. fractal routines.
  5.  
  6. --------------------------------------------------------------------   */
  7.  
  8.  
  9. #include <stdio.h>
  10. #include <stdlib.h>
  11. #include <float.h>
  12. #include <limits.h>
  13. #include <string.h>
  14. #ifdef __TURBOC__
  15. #include <alloc.h>
  16. #elif !defined(__386BSD__)
  17. #include <malloc.h>
  18. #endif
  19. #include "fractint.h"
  20. #include "helpdefs.h"
  21. #include "fractype.h"
  22. #include "prototyp.h"
  23. #include "biginit.h"
  24.  
  25. int bf_math = 0;
  26.  
  27. #if (_MSC_VER >= 700)
  28. #pragma code_seg ("bigsetup_text")     /* place following in an overlay */
  29. #endif
  30.  
  31. /* #define DEBUG */
  32. /* following are handy debugging tools */
  33. #ifdef DEBUG
  34.  
  35. /**********************************************************************/
  36. void show_var_bn(char *s, bn_t n)
  37.     {
  38.         char msg[200];
  39.         strcpy(msg,s);
  40.         strcat(msg," ");
  41.         bntostr(msg+strlen(s),40,n);
  42.         msg[79] = 0;
  43.         stopmsg(0,(char far *)msg);
  44.     }
  45.  
  46. void showcornersdbl(char *s)
  47. {
  48.    char msg[400];
  49.    sprintf(msg,"%s\n\
  50. xxmin= %.20f xxmax= %.20f\n\
  51. yymin= %.20f yymax= %.20f\n\
  52. xx3rd= %.20f yy3rd= %.20f\n\
  53. delxx= %.20Lf delyy= %.20Lf\n\
  54. delx2= %.20Lf dely2= %.20Lf",s,xxmin,xxmax,yymin,yymax,xx3rd,yy3rd,
  55.    delxx, delyy,delxx2, delyy2);
  56.    if(stopmsg(0,msg)==-1)
  57.       goodbye();
  58. }
  59.  
  60. /* show floating point and bignumber corners */
  61. void showcorners(char *s)
  62. {
  63.    int dec=20;
  64.    char msg[100],msg1[100],msg3[100];
  65.    bntostr(msg,dec,bnxmin);
  66.    sprintf(msg1,"bnxmin=%s\nxxmin= %.20f\n\n",msg,xxmin);
  67.    strcpy(msg3,s);
  68.    strcat(msg3,"\n");
  69.    strcat(msg3,msg1);
  70.    bntostr(msg,dec,bnxmax);
  71.    sprintf(msg1,"bnxmax=%s\nxxmax= %.20f\n\n",msg,xxmax);
  72.    strcat(msg3,msg1);
  73.    bntostr(msg,dec,bnymin);
  74.    sprintf(msg1,"bnymin=%s\nyymin= %.20f\n\n",msg,yymin);
  75.    strcat(msg3,msg1);
  76.    bntostr(msg,dec,bnymax);
  77.    sprintf(msg1,"bnymax=%s\nyymax= %.20f\n\n",msg,yymax);
  78.    strcat(msg3,msg1);
  79.    bntostr(msg,dec,bnx3rd);
  80.    sprintf(msg1,"bnx3rd=%s\nxx3rd= %.20f\n\n",msg,xx3rd);
  81.    strcat(msg3,msg1);
  82.    bntostr(msg,dec,bny3rd);
  83.    sprintf(msg1,"bny3rd=%s\nyy3rd= %.20f\n\n",msg,yy3rd);
  84.    strcat(msg3,msg1);
  85.    if(stopmsg(0,msg3)==-1)
  86.       goodbye();
  87. }
  88.  
  89. /* show globals */
  90. void showbfglobals(char *s)
  91. {
  92.    char msg[300];
  93.    sprintf(msg, "%s\n\
  94. bnstep=%d bnlength=%d intlength=%d rlength=%d padding=%d\n\
  95. shiftfactor=%d decimals=%d bflength=%d rbflength=%d \n\
  96. bfshiftfactor=%d bfdecimals=%d ",
  97.                s, bnstep, bnlength, intlength, rlength, padding, 
  98.                shiftfactor, decimals, bflength, rbflength, 
  99.                bfshiftfactor, bfdecimals); 
  100.    if(stopmsg(0,msg)==-1)
  101.       goodbye();
  102. }
  103.  
  104. void showcornersbf(char *s)
  105. {
  106.    int dec=decimals;
  107.    char msg[100],msg1[100],msg3[600];
  108.    if(dec > 20) dec = 20;
  109.    bftostr(msg,dec,bfxmin);
  110.    sprintf(msg1,"bfxmin=%s\nxxmin= %.20f decimals %d bflength %d\n\n",
  111.        msg,xxmin,decimals,bflength);
  112.    strcpy(msg3,s);
  113.    strcat(msg3,"\n");
  114.    strcat(msg3,msg1);
  115.    bftostr(msg,dec,bfxmax);
  116.    sprintf(msg1,"bfxmax=%s\nxxmax= %.20f\n\n",msg,xxmax);
  117.    strcat(msg3,msg1);
  118.    bftostr(msg,dec,bfymin);
  119.    sprintf(msg1,"bfymin=%s\nyymin= %.20f\n\n",msg,yymin);
  120.    strcat(msg3,msg1);
  121.    bftostr(msg,dec,bfymax);
  122.    sprintf(msg1,"bfymax=%s\nyymax= %.20f\n\n",msg,yymax);
  123.    strcat(msg3,msg1);
  124.    bftostr(msg,dec,bfx3rd);
  125.    sprintf(msg1,"bfx3rd=%s\nxx3rd= %.20f\n\n",msg,xx3rd);
  126.    strcat(msg3,msg1);
  127.    bftostr(msg,dec,bfy3rd);
  128.    sprintf(msg1,"bfy3rd=%s\nyy3rd= %.20f\n\n",msg,yy3rd);
  129.    strcat(msg3,msg1);
  130.    if(stopmsg(0,msg3)==-1)
  131.       goodbye();
  132. }
  133.  
  134. void showcornersbfs(char *s)
  135. {
  136.    int dec=20;
  137.    char msg[100],msg1[100],msg3[500];
  138.    bftostr(msg,dec,bfsxmin);
  139.    sprintf(msg1,"bfsxmin=%s\nxxmin= %.20f\n\n",msg,xxmin);
  140.    strcpy(msg3,s);
  141.    strcat(msg3,"\n");
  142.    strcat(msg3,msg1);
  143.    bftostr(msg,dec,bfsxmax);
  144.    sprintf(msg1,"bfsxmax=%s\nxxmax= %.20f\n\n",msg,xxmax);
  145.    strcat(msg3,msg1);
  146.    bftostr(msg,dec,bfsymin);
  147.    sprintf(msg1,"bfsymin=%s\nyymin= %.20f\n\n",msg,yymin);
  148.    strcat(msg3,msg1);
  149.    bftostr(msg,dec,bfsymax);
  150.    sprintf(msg1,"bfsymax=%s\nyymax= %.20f\n\n",msg,yymax);
  151.    strcat(msg3,msg1);
  152.    bftostr(msg,dec,bfsx3rd);
  153.    sprintf(msg1,"bfsx3rd=%s\nxx3rd= %.20f\n\n",msg,xx3rd);
  154.    strcat(msg3,msg1);
  155.    bftostr(msg,dec,bfsy3rd);
  156.    sprintf(msg1,"bfsy3rd=%s\nyy3rd= %.20f\n\n",msg,yy3rd);
  157.    strcat(msg3,msg1);
  158.    if(stopmsg(0,msg3)==-1)
  159.       goodbye();
  160. }
  161.  
  162. void show_two_bf(char *s1,bf_t t1,char *s2, bf_t t2, int digits)
  163. {
  164.    char msg1[200],msg2[200], msg3[400];
  165.    bftostr_e(msg1,digits,t1);
  166.    bftostr_e(msg2,digits,t2);
  167.    sprintf(msg3,"\n%s->%s\n%s->%s",s1,msg1,s2,msg2);
  168.    if(stopmsg(0,msg3)==-1)
  169.       goodbye();
  170. }
  171.  
  172. void show_three_bf(char *s1,bf_t t1,char *s2, bf_t t2, char *s3, bf_t t3, int digits)
  173. {
  174.    char msg1[200],msg2[200], msg3[200], msg4[600];
  175.    bftostr_e(msg1,digits,t1);
  176.    bftostr_e(msg2,digits,t2);
  177.    bftostr_e(msg3,digits,t3);
  178.    sprintf(msg4,"\n%s->%s\n%s->%s\n%s->%s",s1,msg1,s2,msg2,s3,msg3);
  179.    if(stopmsg(0,msg4)==-1)
  180.       goodbye();
  181. }
  182.  
  183. /* for aspect ratio debugging */
  184. void showaspect(char *s)
  185. {
  186.    bf_t bt1,bt2,aspect;
  187.    char msg[100],str[100];
  188.    int saved; saved = save_stack();
  189.    bt1    = alloc_stack(rbflength+2);
  190.    bt2    = alloc_stack(rbflength+2);
  191.    aspect = alloc_stack(rbflength+2);
  192.    sub_bf(bt1,bfxmax,bfxmin);
  193.    sub_bf(bt2,bfymax,bfymin);
  194.    div_bf(aspect,bt2,bt1);
  195.    bftostr(str,10,aspect); 
  196.    sprintf(msg,"aspect %s\nfloat %13.10f\nbf    %s\n\n",
  197.             s,
  198.             (yymax-yymin)/(xxmax-xxmin),
  199.             str);
  200.    if(stopmsg(0,msg)==-1)
  201.       goodbye();
  202.    restore_stack(saved);
  203. }
  204.  
  205. /* compare a double and bignumber */
  206. void comparevalues(char *s, LDBL x, bn_t bnx)
  207. {
  208.    int dec=40;
  209.    char msg[100],msg1[100];
  210.    bntostr(msg,dec,bnx);
  211.    sprintf(msg1,"%s\nbignum=%s\ndouble=%.20Lf\n\n",s,msg,x);
  212.    if(stopmsg(0,msg1)==-1)
  213.       goodbye();
  214. }
  215. /* compare a double and bignumber */
  216. void comparevaluesbf(char *s, LDBL x, bf_t bfx)
  217. {
  218.    int dec=40;
  219.    char msg[300],msg1[300];
  220.    bftostr_e(msg,dec,bfx);
  221.    sprintf(msg1,"%s\nbignum=%s\ndouble=%.20Lf\n\n",s,msg,x);
  222.    if(stopmsg(0,msg1)==-1)
  223.       goodbye();
  224. }
  225.  
  226. /**********************************************************************/
  227. void show_var_bf(char *s, bf_t n)
  228.     {
  229.         char msg[200];
  230.         strcpy(msg,s);
  231.         strcat(msg," ");
  232.         bftostr_e(msg+strlen(s),40,n);
  233.         msg[79] = 0;
  234.         if(stopmsg(0,msg)==-1)
  235.             goodbye();
  236.     }
  237.  
  238. #endif
  239.  
  240. void bfcornerstofloat(void)
  241. {
  242.    int i;
  243.    if(bf_math)
  244.    {
  245.       xxmin = (double)bftofloat(bfxmin);
  246.       yymin = (double)bftofloat(bfymin);
  247.       xxmax = (double)bftofloat(bfxmax);
  248.       yymax = (double)bftofloat(bfymax);
  249.       xx3rd = (double)bftofloat(bfx3rd);
  250.       yy3rd = (double)bftofloat(bfy3rd);
  251.    }
  252.    for(i=0;i<MAXPARAMS;i++)
  253.       if(typehasparm(fractype,i))
  254.          param[i] = (double)bftofloat(bfparms[i]);
  255. }
  256.  
  257. #if (_MSC_VER >= 700)
  258. #pragma code_seg ( )     /* back to normal segment */
  259. #endif
  260.  
  261. /* -------------------------------------------------------------------- */
  262. /*    Bignumber Bailout Routines                                        */
  263. /* -------------------------------------------------------------------- */
  264.  
  265. /* mandel_bntoint() can only be used for intlength of 1 */
  266. #define mandel_bntoint(n) (*(n + bnlength - 1)) /* assumes intlength of 1 */
  267.  
  268. /* Note:                                             */
  269. /* No need to set magnitude                          */
  270. /* as color schemes that need it calculate it later. */
  271.  
  272. int near bnMODbailout()
  273. {
  274.    long longmagnitude;
  275.  
  276.    square_bn(bntmpsqrx, bnnew.x);
  277.    square_bn(bntmpsqry, bnnew.y);
  278.    add_bn(bntmp, bntmpsqrx+shiftfactor, bntmpsqry+shiftfactor);
  279.  
  280.    longmagnitude = bntoint(bntmp);  /* works with any fractal type */
  281.    if (longmagnitude >= (long)rqlim)
  282.       return 1;
  283.    copy_bn(bnold.x, bnnew.x);
  284.    copy_bn(bnold.y, bnnew.y);
  285.    return 0;
  286. }
  287.   
  288. int near bnREALbailout()
  289. {
  290.    long longtempsqrx;
  291.  
  292.    square_bn(bntmpsqrx, bnnew.x);
  293.    square_bn(bntmpsqry, bnnew.y);
  294.    longtempsqrx = bntoint(bntmpsqrx+shiftfactor);
  295.    if (longtempsqrx >= (long)rqlim)
  296.       return 1;
  297.    copy_bn(bnold.x, bnnew.x);
  298.    copy_bn(bnold.y, bnnew.y);
  299.    return 0;
  300. }
  301.  
  302.  
  303. int near bnIMAGbailout()
  304. {
  305.    long longtempsqry;
  306.  
  307.    square_bn(bntmpsqrx, bnnew.x);
  308.    square_bn(bntmpsqry, bnnew.y);
  309.    longtempsqry = bntoint(bntmpsqry+shiftfactor);
  310.    if (longtempsqry >= (long)rqlim)
  311.       return 1;
  312.    copy_bn(bnold.x, bnnew.x);
  313.    copy_bn(bnold.y, bnnew.y);
  314.    return(0);
  315. }
  316.  
  317. int near bnORbailout()
  318. {
  319.    long longtempsqrx, longtempsqry;
  320.  
  321.    square_bn(bntmpsqrx, bnnew.x);
  322.    square_bn(bntmpsqry, bnnew.y);
  323.    longtempsqrx = bntoint(bntmpsqrx+shiftfactor);
  324.    longtempsqry = bntoint(bntmpsqry+shiftfactor);
  325.    if(longtempsqrx >= (long)rqlim || longtempsqry >= (long)rqlim)
  326.       return 1;
  327.    copy_bn(bnold.x, bnnew.x);
  328.    copy_bn(bnold.y, bnnew.y);
  329.    return(0);
  330. }
  331.  
  332. int near bnANDbailout()
  333. {
  334.    long longtempsqrx, longtempsqry;
  335.  
  336.    square_bn(bntmpsqrx, bnnew.x);
  337.    square_bn(bntmpsqry, bnnew.y);
  338.    longtempsqrx = bntoint(bntmpsqrx+shiftfactor); 
  339.    longtempsqry = bntoint(bntmpsqry+shiftfactor); 
  340.    if(longtempsqrx >= (long)rqlim && longtempsqry >= (long)rqlim)
  341.       return 1;
  342.    copy_bn(bnold.x, bnnew.x);
  343.    copy_bn(bnold.y, bnnew.y);
  344.    return(0);
  345. }
  346.  
  347. int near bfMODbailout()
  348. {
  349.    long longmagnitude;
  350.  
  351.    square_bf(bftmpsqrx, bfnew.x);
  352.    square_bf(bftmpsqry, bfnew.y);
  353.    add_bf(bftmp, bftmpsqrx, bftmpsqry);
  354.  
  355.    longmagnitude = bftoint(bftmp);
  356.    if (longmagnitude >= (long)rqlim)
  357.       return 1;
  358.    copy_bf(bfold.x, bfnew.x);
  359.    copy_bf(bfold.y, bfnew.y);
  360.    return 0;
  361. }
  362.   
  363. int near bfREALbailout()
  364. {
  365.    long longtempsqrx;
  366.  
  367.    square_bf(bftmpsqrx, bfnew.x);
  368.    square_bf(bftmpsqry, bfnew.y);
  369.    longtempsqrx = bftoint(bftmpsqrx);
  370.    if (longtempsqrx >= (long)rqlim)
  371.       return 1;
  372.    copy_bf(bfold.x, bfnew.x);
  373.    copy_bf(bfold.y, bfnew.y);
  374.    return 0;
  375. }
  376.  
  377.  
  378. int near bfIMAGbailout()
  379. {
  380.    long longtempsqry;
  381.  
  382.    square_bf(bftmpsqrx, bfnew.x);
  383.    square_bf(bftmpsqry, bfnew.y);
  384.    longtempsqry = bftoint(bftmpsqry);
  385.    if (longtempsqry >= (long)rqlim)
  386.       return 1;
  387.    copy_bf(bfold.x, bfnew.x);
  388.    copy_bf(bfold.y, bfnew.y);
  389.    return(0);
  390. }
  391.  
  392. int near bfORbailout()
  393. {
  394.    long longtempsqrx, longtempsqry;
  395.  
  396.    square_bf(bftmpsqrx, bfnew.x);
  397.    square_bf(bftmpsqry, bfnew.y);
  398.    longtempsqrx = bftoint(bftmpsqrx);
  399.    longtempsqry = bftoint(bftmpsqry);
  400.    if(longtempsqrx >= (long)rqlim || longtempsqry >= (long)rqlim)
  401.       return 1;
  402.    copy_bf(bfold.x, bfnew.x);
  403.    copy_bf(bfold.y, bfnew.y);
  404.    return(0);
  405. }
  406.  
  407. int near bfANDbailout()
  408. {
  409.    long longtempsqrx, longtempsqry;
  410.  
  411.    square_bf(bftmpsqrx, bfnew.x);
  412.    square_bf(bftmpsqry, bfnew.y);
  413.    longtempsqrx = bftoint(bftmpsqrx);
  414.    longtempsqry = bftoint(bftmpsqry);
  415.    if(longtempsqrx >= (long)rqlim && longtempsqry >= (long)rqlim)
  416.       return 1;
  417.    copy_bf(bfold.x, bfnew.x);
  418.    copy_bf(bfold.y, bfnew.y);
  419.    return(0);
  420. }
  421.  
  422. #if (_MSC_VER >= 700)
  423. #pragma code_seg ("bigsetup_text")     /* place following in an overlay */
  424. #endif
  425.  
  426. int MandelbnSetup()
  427. {
  428.    /* this should be set up dynamically based on corners */
  429.    bn_t bntemp1, bntemp2;
  430.    int saved; saved = save_stack();
  431.    bntemp1 = alloc_stack(bnlength);
  432.    bntemp2 = alloc_stack(bnlength);
  433.  
  434.    bftobn(bnxmin,bfxmin);
  435.    bftobn(bnxmax,bfxmax);
  436.    bftobn(bnymin,bfymin);
  437.    bftobn(bnymax,bfymax);
  438.    bftobn(bnx3rd,bfx3rd);
  439.    bftobn(bny3rd,bfy3rd);
  440.  
  441.    bf_math = BIGNUM;
  442.  
  443.    /* bnxdel = (bnxmax - bnx3rd)/(xdots-1) */
  444.    sub_bn(bnxdel, bnxmax, bnx3rd);
  445.    div_a_bn_int(bnxdel, (U16)(xdots - 1));
  446.  
  447.    /* bnydel = (bnymax - bny3rd)/(ydots-1) */
  448.    sub_bn(bnydel, bnymax, bny3rd);
  449.    div_a_bn_int(bnydel, (U16)(ydots - 1));
  450.  
  451.    /* bnxdel2 = (bnx3rd - bnxmin)/(ydots-1) */
  452.    sub_bn(bnxdel2, bnx3rd, bnxmin);
  453.    div_a_bn_int(bnxdel2, (U16)(ydots - 1));
  454.  
  455.    /* bnydel2 = (bny3rd - bnymin)/(xdots-1) */
  456.    sub_bn(bnydel2, bny3rd, bnymin);
  457.    div_a_bn_int(bnydel2, (U16)(xdots - 1));
  458.  
  459.    abs_bn(bnclosenuff,bnxdel);
  460.    if(cmp_bn(abs_bn(bntemp1,bnxdel2),bnclosenuff) > 0)
  461.       copy_bn(bnclosenuff,bntemp1);
  462.    if(cmp_bn(abs_bn(bntemp1,bnydel),abs_bn(bntemp2,bnydel2)) > 0)
  463.    {
  464.       if(cmp_bn(bntemp1,bnclosenuff) > 0)
  465.          copy_bn(bnclosenuff,bntemp1);
  466.    } 
  467.    else if(cmp_bn(bntemp2,bnclosenuff) > 0)
  468.       copy_bn(bnclosenuff,bntemp2);
  469.    {
  470.       int t;
  471.       t = abs(periodicitycheck);
  472.       while(t--)
  473.          half_a_bn(bnclosenuff);
  474.    }
  475.  
  476.    c_exp = (int)param[2];
  477.    switch (fractype)
  478.    {
  479.       case JULIAFP:
  480.          bftobn(bnparm.x, bfparms[0]);
  481.          bftobn(bnparm.y, bfparms[1]);
  482.          break;
  483.       case FPMANDELZPOWER:
  484.          init_big_pi();
  485.          if((double)c_exp == param[2] && (c_exp & 1)) /* odd exponents */
  486.             symmetry = XYAXIS_NOPARM;
  487.          if(param[3] != 0)
  488.             symmetry = NOSYM;
  489.          break;
  490.       case FPJULIAZPOWER:
  491.          init_big_pi();
  492.          bftobn(bnparm.x, bfparms[0]);
  493.          bftobn(bnparm.y, bfparms[1]);
  494.          if((c_exp & 1) || param[3] != 0.0 || (double)c_exp != param[2] )
  495.             symmetry = NOSYM;
  496.          break;
  497.    }
  498.  
  499. /* at the present time, parameters are kept in float, but want to keep
  500.       the arbitrary precision logic intact. The next two lines, if used,
  501.       would disguise and breaking of the arbitrary precision logic */
  502.    /*
  503.    floattobn(bnparm.x,param[0]);
  504.    floattobn(bnparm.y,param[1]);
  505.    */
  506.    restore_stack(saved);
  507.    return (1);
  508. }
  509.  
  510. int MandelbfSetup()
  511. {
  512.    /* this should be set up dynamically based on corners */
  513.    bf_t bftemp1, bftemp2;
  514.    int saved; saved = save_stack();
  515.    bftemp1 = alloc_stack(bflength+2);
  516.    bftemp2 = alloc_stack(bflength+2);
  517.  
  518.    bf_math = BIGFLT;
  519.  
  520.    /* bfxdel = (bfxmax - bfx3rd)/(xdots-1) */
  521.    sub_bf(bfxdel, bfxmax, bfx3rd);
  522.    div_a_bf_int(bfxdel, (U16)(xdots - 1));
  523.  
  524.    /* bfydel = (bfymax - bfy3rd)/(ydots-1) */
  525.    sub_bf(bfydel, bfymax, bfy3rd);
  526.    div_a_bf_int(bfydel, (U16)(ydots - 1));
  527.  
  528.    /* bfxdel2 = (bfx3rd - bfxmin)/(ydots-1) */
  529.    sub_bf(bfxdel2, bfx3rd, bfxmin);
  530.    div_a_bf_int(bfxdel2, (U16)(ydots - 1));
  531.  
  532.    /* bfydel2 = (bfy3rd - bfymin)/(xdots-1) */
  533.    sub_bf(bfydel2, bfy3rd, bfymin);
  534.    div_a_bf_int(bfydel2, (U16)(xdots - 1));
  535.  
  536.    abs_bf(bfclosenuff,bfxdel);
  537.    if(cmp_bf(abs_bf(bftemp1,bfxdel2),bfclosenuff) > 0)
  538.       copy_bf(bfclosenuff,bftemp1);
  539.    if(cmp_bf(abs_bf(bftemp1,bfydel),abs_bf(bftemp2,bfydel2)) > 0)
  540.    {
  541.       if(cmp_bf(bftemp1,bfclosenuff) > 0)
  542.          copy_bf(bfclosenuff,bftemp1);
  543.    }
  544.    else if(cmp_bf(bftemp2,bfclosenuff) > 0)
  545.       copy_bf(bfclosenuff,bftemp2);
  546.    {
  547.       int t;
  548.       t = abs(periodicitycheck);
  549.       while(t--)
  550.          half_a_bf(bfclosenuff);
  551.    }
  552.  
  553.    c_exp = (int)param[2];
  554.    switch (fractype)
  555.    {
  556.       case JULIAFP:
  557.          copy_bf(bfparm.x, bfparms[0]);
  558.          copy_bf(bfparm.y, bfparms[1]);
  559.          break;
  560.       case FPMANDELZPOWER:
  561.          init_big_pi();
  562.          if((double)c_exp == param[2] && (c_exp & 1)) /* odd exponents */
  563.             symmetry = XYAXIS_NOPARM;
  564.          if(param[3] != 0)
  565.             symmetry = NOSYM;
  566.          break;
  567.       case FPJULIAZPOWER:
  568.          init_big_pi();
  569.          copy_bf(bfparm.x, bfparms[0]);
  570.          copy_bf(bfparm.y, bfparms[1]);
  571.          if((c_exp & 1) || param[3] != 0.0 || (double)c_exp != param[2] )
  572.             symmetry = NOSYM;
  573.          break;
  574.    }
  575.  
  576.    restore_stack(saved);
  577.    return (1);
  578. }
  579.  
  580. #if (_MSC_VER >= 700)
  581. #pragma code_seg ( )     /* back to normal segment */
  582. #endif
  583.  
  584. int mandelbn_per_pixel()
  585. {
  586.    /* parm.x = xxmin + col*delx + row*delx2 */
  587.    mult_bn_int(bnparm.x, bnxdel, (U16)col);
  588.    mult_bn_int(bntmp, bnxdel2, (U16)row);
  589.  
  590.    add_a_bn(bnparm.x, bntmp);
  591.    add_a_bn(bnparm.x, bnxmin);
  592.  
  593.    /* parm.y = yymax - row*dely - col*dely2; */
  594.    /* note: in next four lines, bnold is just used as a temporary variable */
  595.    mult_bn_int(bnold.x, bnydel,  (U16)row);
  596.    mult_bn_int(bnold.y, bnydel2, (U16)col);
  597.    add_a_bn(bnold.x, bnold.y);
  598.    sub_bn(bnparm.y, bnymax, bnold.x);
  599.  
  600.    copy_bn(bnold.x, bnparm.x);
  601.    copy_bn(bnold.y, bnparm.y);
  602.  
  603.    if(inside == -60 || inside == -61)
  604.    {
  605.       /* kludge to match "Beauty of Fractals" picture since we start
  606.      Mandelbrot iteration with init rather than 0 */
  607.       floattobn(bnold.x,param[0]); /* initial pertubation of parameters set */
  608.       floattobn(bnold.y,param[1]);
  609.       coloriter = -1;
  610.    }
  611.    else
  612.    {
  613.      floattobn(bnnew.x,param[0]); 
  614.      floattobn(bnnew.y,param[1]);
  615.      add_a_bn(bnold.x,bnnew.x);
  616.      add_a_bn(bnold.y,bnnew.y);
  617.    }
  618.  
  619.    /* square has side effect - must copy first */
  620.    copy_bn(bnnew.x, bnold.x);
  621.    copy_bn(bnnew.y, bnold.y);
  622.  
  623.    /* Square these to rlength bytes of precision */
  624.    square_bn(bntmpsqrx, bnnew.x);
  625.    square_bn(bntmpsqry, bnnew.y);
  626.  
  627.    return (1);            /* 1st iteration has been done */
  628. }
  629.  
  630. int mandelbf_per_pixel()
  631. {
  632.    /* parm.x = xxmin + col*delx + row*delx2 */
  633.    mult_bf_int(bfparm.x, bfxdel, (U16)col);
  634.    mult_bf_int(bftmp, bfxdel2, (U16)row);
  635.  
  636.    add_a_bf(bfparm.x, bftmp);
  637.    add_a_bf(bfparm.x, bfxmin);
  638.  
  639.    /* parm.y = yymax - row*dely - col*dely2; */
  640.    /* note: in next four lines, bfold is just used as a temporary variable */
  641.    mult_bf_int(bfold.x, bfydel,  (U16)row);
  642.    mult_bf_int(bfold.y, bfydel2, (U16)col);
  643.    add_a_bf(bfold.x, bfold.y);
  644.    sub_bf(bfparm.y, bfymax, bfold.x);
  645.  
  646.    copy_bf(bfold.x, bfparm.x);
  647.    copy_bf(bfold.y, bfparm.y);
  648.  
  649.    if(inside == -60 || inside == -61)
  650.    {
  651.       /* kludge to match "Beauty of Fractals" picture since we start
  652.      Mandelbrot iteration with init rather than 0 */
  653.       floattobf(bfold.x,param[0]); /* initial pertubation of parameters set */
  654.       floattobf(bfold.y,param[1]);
  655.       coloriter = -1;
  656.    }
  657.    else
  658.    {
  659.      floattobf(bfnew.x,param[0]);
  660.      floattobf(bfnew.y,param[1]);
  661.      add_a_bf(bfold.x,bfnew.x);
  662.      add_a_bf(bfold.y,bfnew.y);
  663.    }
  664.  
  665.    /* square has side effect - must copy first */
  666.    copy_bf(bfnew.x, bfold.x);
  667.    copy_bf(bfnew.y, bfold.y);
  668.  
  669.    /* Square these to rbflength bytes of precision */
  670.    square_bf(bftmpsqrx, bfnew.x);
  671.    square_bf(bftmpsqry, bfnew.y);
  672.  
  673.    return (1);            /* 1st iteration has been done */
  674. }
  675.  
  676. juliabn_per_pixel()
  677. {
  678.    /* old.x = xxmin + col*delx + row*delx2 */
  679.    mult_bn_int(bnold.x, bnxdel, (U16)col);
  680.    mult_bn_int(bntmp, bnxdel2, (U16)row);
  681.  
  682.    add_a_bn(bnold.x, bntmp);
  683.    add_a_bn(bnold.x, bnxmin);
  684.  
  685.    /* old.y = yymax - row*dely - col*dely2; */
  686.    /* note: in next four lines, bnnew is just used as a temporary variable */
  687.    mult_bn_int(bnnew.x, bnydel,  (U16)row);
  688.    mult_bn_int(bnnew.y, bnydel2, (U16)col);
  689.    add_a_bn(bnnew.x, bnnew.y);
  690.    sub_bn(bnold.y, bnymax, bnnew.x);
  691.    
  692.    /* square has side effect - must copy first */
  693.    copy_bn(bnnew.x, bnold.x);
  694.    copy_bn(bnnew.y, bnold.y);
  695.  
  696.    /* Square these to rlength bytes of precision */
  697.    square_bn(bntmpsqrx, bnnew.x);
  698.    square_bn(bntmpsqry, bnnew.y);
  699.  
  700.    return (1);            /* 1st iteration has been done */
  701. }
  702.  
  703. juliabf_per_pixel()
  704. {
  705.    /* old.x = xxmin + col*delx + row*delx2 */
  706.    mult_bf_int(bfold.x, bfxdel, (U16)col);
  707.    mult_bf_int(bftmp, bfxdel2, (U16)row);
  708.  
  709.    add_a_bf(bfold.x, bftmp);
  710.    add_a_bf(bfold.x, bfxmin);
  711.  
  712.    /* old.y = yymax - row*dely - col*dely2; */
  713.    /* note: in next four lines, bfnew is just used as a temporary variable */
  714.    mult_bf_int(bfnew.x, bfydel,  (U16)row);
  715.    mult_bf_int(bfnew.y, bfydel2, (U16)col);
  716.    add_a_bf(bfnew.x, bfnew.y);
  717.    sub_bf(bfold.y, bfymax, bfnew.x);
  718.    
  719.    /* square has side effect - must copy first */
  720.    copy_bf(bfnew.x, bfold.x);
  721.    copy_bf(bfnew.y, bfold.y);
  722.  
  723.    /* Square these to rbflength bytes of precision */
  724.    square_bf(bftmpsqrx, bfnew.x);
  725.    square_bf(bftmpsqry, bfnew.y);
  726.  
  727.    return (1);            /* 1st iteration has been done */
  728. }
  729.  
  730.  
  731. JuliabnFractal()
  732. {
  733.    /* Don't forget, with bn_t numbers, after multiplying or squaring */
  734.    /* you must shift over by shiftfactor to get the bn number.          */
  735.  
  736.    /* bntmpsqrx and bntmpsqry were previously squared before getting to */
  737.    /* this function, so they must be shifted.                           */
  738.  
  739.    /* new.x = tmpsqrx - tmpsqry + parm.x;   */
  740.    sub_a_bn(bntmpsqrx+shiftfactor, bntmpsqry+shiftfactor);
  741.    add_bn(bnnew.x, bntmpsqrx+shiftfactor, bnparm.x);
  742.  
  743.    /* new.y = 2 * bnold.x * bnold.y + parm.y; */
  744.    mult_bn(bntmp, bnold.x, bnold.y); /* ok to use unsafe here */
  745.    double_a_bn(bntmp+shiftfactor);
  746.    add_bn(bnnew.y, bntmp+shiftfactor, bnparm.y);
  747.  
  748.    return bignumbailout();
  749. }
  750.  
  751. JuliabfFractal()
  752. {
  753.    /* new.x = tmpsqrx - tmpsqry + parm.x;   */
  754.    sub_a_bf(bftmpsqrx, bftmpsqry);
  755.    add_bf(bfnew.x, bftmpsqrx, bfparm.x);
  756.  
  757.    /* new.y = 2 * bfold.x * bfold.y + parm.y; */
  758.    mult_bf(bftmp, bfold.x, bfold.y); /* ok to use unsafe here */
  759.    double_a_bf(bftmp);
  760.    add_bf(bfnew.y, bftmp, bfparm.y);
  761.    return bigfltbailout();
  762. }
  763.  
  764. JuliaZpowerbnFractal()
  765. {
  766.    _BNCMPLX parm2;
  767.    int saved; saved = save_stack();
  768.  
  769.    parm2.x = alloc_stack(bnlength);
  770.    parm2.y = alloc_stack(bnlength);
  771.  
  772.    floattobn(parm2.x,param[2]);
  773.    floattobn(parm2.y,param[3]);
  774.    ComplexPower_bn(&bnnew,&bnold,&parm2); 
  775.    add_bn(bnnew.x,bnparm.x,bnnew.x+shiftfactor);
  776.    add_bn(bnnew.y,bnparm.y,bnnew.y+shiftfactor);
  777.    restore_stack(saved);
  778.    return bignumbailout();
  779. }
  780.  
  781. JuliaZpowerbfFractal()
  782. {
  783.    _BFCMPLX parm2;
  784.    int saved; saved = save_stack();
  785.  
  786.    parm2.x = alloc_stack(bflength+2);
  787.    parm2.y = alloc_stack(bflength+2);
  788.  
  789.    floattobf(parm2.x,param[2]);
  790.    floattobf(parm2.y,param[3]);
  791.    ComplexPower_bf(&bfnew,&bfold,&parm2);
  792.    add_bf(bfnew.x,bfparm.x,bfnew.x);
  793.    add_bf(bfnew.y,bfparm.y,bfnew.y);
  794.    restore_stack(saved);
  795.    return bigfltbailout();
  796. }
  797.  
  798.  
  799. #if 0
  800. /*
  801. the following is an example of how you can take advantage of the bn_t
  802. format to squeeze a little more precision out of the calculations.
  803. */
  804. JuliabnFractal()
  805. {
  806.    int oldbnlength;
  807.    bn_t mod;
  808.    /* using partial precision multiplications */
  809.  
  810.    /* bnnew.x = bntmpsqrx - bntmpsqry + bnparm.x;   */
  811.    /*
  812.     * Since tmpsqrx and tmpsqry where just calculated to rlength bytes of
  813.     * precision, we might as well keep that extra precision in this next
  814.     * subtraction.  Therefore, use rlength as the length.
  815.     */
  816.     
  817.    oldbnlength = bnlength;
  818.    bnlength = rlength; sub_a_bn(bntmpsqrx, bntmpsqry); bnlength = oldbnlength;
  819.  
  820.    /*
  821.     * Now that bntmpsqry has been sutracted from bntmpsqrx, we need to treat
  822.     * tmpsqrx as a single width bignumber, so shift to bntmpsqrx+shiftfactor.
  823.     */
  824.    add_bn(bnnew.x, bntmpsqrx + shiftfactor, bnparm.x);
  825.  
  826.    /* new.y = 2 * bnold.x * bnold.y + old.y; */
  827.    /* Multiply bnold.x*bnold.y to rlength precision. */
  828.    mult_bn(bntmp, bnold.x, bnold.y);
  829.  
  830.    /*
  831.     * Double bnold.x*bnold.y by shifting bits, including one of those bits
  832.     * calculated in the previous mult_bn().  Therefore, use rlength.
  833.     */
  834.    bnlength = rlength; double_a_bn(bntmp); bnlength = oldbnlength;
  835.  
  836.    /* Convert back to a single width bignumber and add bnparm.y */
  837.    add_bn(bnnew.y, bntmp + shiftfactor, bnparm.y);
  838.  
  839.    copy_bn(bnold.x, bnnew.x);
  840.    copy_bn(bnold.y, bnnew.y);
  841.  
  842.    /* Square these to rlength bytes of precision */
  843.    square_bn(bntmpsqrx, bnold.x);
  844.    square_bn(bntmpsqry, bnold.y);
  845.  
  846.    /* And add the full rlength precision to get those extra bytes */
  847.    bnlength = rlength;
  848.    add_bn(bntmp, bntmpsqrx, bntmpsqry);
  849.    bnlength = oldbnlength;
  850.  
  851.    mod = bntmp + (rlength) - (intlength << 1);    /* where int part starts
  852.                          * after mult */
  853.    /*
  854.     * equivalent to, but faster than, mod = bn_int(tmp+shiftfactor);
  855.     */
  856.    
  857.    if ((magnitude = *mod) >= rqlim)
  858.       return (1);
  859.    return (0);
  860. }
  861. #endif
  862.  
  863. _CMPLX cmplxbntofloat(_BNCMPLX *s)
  864. {
  865.    _CMPLX t;
  866.    t.x = (double)bntofloat(s->x);
  867.    t.y = (double)bntofloat(s->y);
  868.    return(t);
  869. }
  870.  
  871. _CMPLX cmplxbftofloat(_BFCMPLX *s)
  872. {
  873.    _CMPLX t; 
  874.    t.x = (double)bftofloat(s->x);
  875.    t.y = (double)bftofloat(s->y);
  876.    return(t);
  877. }
  878.  
  879. _BFCMPLX *cmplxlog_bf(_BFCMPLX *t, _BFCMPLX *s)
  880. {
  881.    square_bf(t->x,s->x); 
  882.    square_bf(t->y,s->y);
  883.    add_a_bf(t->x,t->y);
  884.    ln_bf(t->x,t->x);
  885.    half_a_bf(t->x);
  886.    atan2_bf(t->y,s->y,s->x);
  887.    return(t);
  888. }
  889.  
  890. _BFCMPLX *cplxmul_bf( _BFCMPLX *t, _BFCMPLX *x, _BFCMPLX *y)
  891. {
  892.    bf_t tmp1;
  893.    int saved; saved = save_stack();
  894.    tmp1 = alloc_stack(rbflength+2);
  895.    mult_bf(t->x, x->x, y->x);
  896.    mult_bf(t->y, x->y, y->y);
  897.    sub_bf(t->x,t->x,t->y); 
  898.  
  899.    mult_bf(tmp1, x->x, y->y);
  900.    mult_bf(t->y, x->y, y->x);
  901.    add_bf(t->y,tmp1,t->y);
  902.    restore_stack(saved);
  903.    return(t);
  904. }
  905.  
  906. _BFCMPLX *ComplexPower_bf(_BFCMPLX *t, _BFCMPLX *xx, _BFCMPLX *yy) 
  907. {
  908.    _BFCMPLX tmp;
  909.    bf_t e2x, siny, cosy;
  910.    int saved; saved = save_stack();
  911.    e2x  = alloc_stack(rbflength+2);
  912.    siny = alloc_stack(rbflength+2);
  913.    cosy = alloc_stack(rbflength+2);
  914.    tmp.x = alloc_stack(rbflength+2);
  915.    tmp.y = alloc_stack(rbflength+2);
  916.  
  917.    /* 0 raised to anything is 0 */
  918.    if (is_bf_zero(xx->x) && is_bf_zero(xx->y))
  919.       {
  920.       clear_bf(t->x);
  921.       clear_bf(t->y);
  922.       return(t);
  923.       }
  924.  
  925.    cmplxlog_bf(t, xx);
  926.    cplxmul_bf(&tmp, t, yy);
  927.    exp_bf(e2x,tmp.x);
  928.    sincos_bf(siny,cosy,tmp.y);
  929.    mult_bf(t->x, e2x, cosy);
  930.    mult_bf(t->y, e2x, siny);
  931.    restore_stack(saved);
  932.    return(t);
  933. }
  934.  
  935. _BNCMPLX *cmplxlog_bn(_BNCMPLX *t, _BNCMPLX *s)
  936. {
  937.    square_bn(t->x,s->x);
  938.    square_bn(t->y,s->y);
  939.    add_a_bn(t->x+shiftfactor,t->y+shiftfactor);
  940.    ln_bn(t->x,t->x+shiftfactor);
  941.    half_a_bn(t->x);
  942.    atan2_bn(t->y,s->y,s->x);
  943.    return(t);
  944. }
  945.  
  946. _BNCMPLX *cplxmul_bn( _BNCMPLX *t, _BNCMPLX *x, _BNCMPLX *y)
  947. {
  948.    bn_t tmp1;
  949.    int saved; saved = save_stack();
  950.    tmp1 = alloc_stack(rlength);
  951.    mult_bn(t->x, x->x, y->x);
  952.    mult_bn(t->y, x->y, y->y);
  953.    sub_bn(t->x,t->x+shiftfactor,t->y+shiftfactor);
  954.  
  955.    mult_bn(tmp1, x->x, y->y);
  956.    mult_bn(t->y, x->y, y->x);
  957.    add_bn(t->y,tmp1+shiftfactor,t->y+shiftfactor);
  958.    restore_stack(saved);
  959.    return(t);
  960. }
  961.  
  962. /* note: ComplexPower_bn() returns need to be +shiftfactor'ed */
  963. _BNCMPLX *ComplexPower_bn(_BNCMPLX *t, _BNCMPLX *xx, _BNCMPLX *yy) 
  964. {
  965.    _BNCMPLX tmp;
  966.    bn_t e2x, siny, cosy;
  967.    int saved; saved = save_stack();
  968.    e2x  = alloc_stack(bnlength);
  969.    siny = alloc_stack(bnlength);
  970.    cosy = alloc_stack(bnlength);
  971.    tmp.x = alloc_stack(rlength);
  972.    tmp.y = alloc_stack(rlength);
  973.  
  974.    /* 0 raised to anything is 0 */
  975.    if (is_bn_zero(xx->x) && is_bn_zero(xx->y))
  976.       {
  977.       clear_bn(t->x);
  978.       clear_bn(t->y);
  979.       return(t);
  980.       }
  981.  
  982.    cmplxlog_bn(t, xx);
  983.    cplxmul_bn(&tmp, t, yy);
  984.    exp_bn(e2x,tmp.x);
  985.    sincos_bn(siny,cosy,tmp.y);
  986.    mult_bn(t->x, e2x, cosy);
  987.    mult_bn(t->y, e2x, siny);
  988.    restore_stack(saved);
  989.    return(t);
  990. }
  991.